EBImage is a Bioconductor package that offers tools to process and analyze images.

Installation

To install EBImage, start R (version “3.5”) and enter:

#if (!requireNamespace("BiocManager", quietly = TRUE))
#    install.packages("BiocManager")
#BiocManager::install("EBImage", version = "3.8")

Loading

To load EBImage, use the following command:

library("EBImage")

Reading

Images are loaded using the function readImage(), which takes as input a character vector of file names or URLs. The loaded images are stored into an Image object.

#files=list.files(system.file("images",package="EBImage"))
#sampleColor=readImage(system.file("images","sample-color.png",package="EBImage"))
#sample=readImage(system.file("images","sample.png",package="EBImage"))
#shapes=readImage(system.file("images","shapes.png",package="EBImage"))
vesuvio=readImage("vesuvio.png")
cells=readImage(system.file("images","cells.tif",package="EBImage"))
nuclei=readImage(system.file("images","nuclei.tif",package="EBImage"))
greyVesuvio=channel(vesuvio,'grey')

Grayscale images

str(greyVesuvio)
## Formal class 'Image' [package "EBImage"] with 2 slots
##   ..@ .Data    : num [1:851, 1:315] 0.639 0.642 0.648 0.648 0.641 ...
##   ..@ colormode: int 0
print(greyVesuvio)
## Image 
##   colorMode    : Grayscale 
##   storage.mode : double 
##   dim          : 851 315 
##   frames.total : 1 
##   frames.render: 1 
## 
## imageData(object)[1:5,1:6]
##           [,1]      [,2]      [,3]      [,4]      [,5]      [,6]
## [1,] 0.6392157 0.6339869 0.6418301 0.6483660 0.6418301 0.6392157
## [2,] 0.6418301 0.6418301 0.6431373 0.6457516 0.6379085 0.6379085
## [3,] 0.6483660 0.6522876 0.6418301 0.6392157 0.6496732 0.6562092
## [4,] 0.6483660 0.6418301 0.6392157 0.6483660 0.6509804 0.6457516
## [5,] 0.6405229 0.6339869 0.6444444 0.6483660 0.6418301 0.6339869
dim(greyVesuvio)
## [1] 851 315
imageData(greyVesuvio)[10:14,4:8]
##           [,1]      [,2]      [,3]      [,4]      [,5]
## [1,] 0.6483660 0.6444444 0.6366013 0.6444444 0.6483660
## [2,] 0.6444444 0.6366013 0.6405229 0.6444444 0.6405229
## [3,] 0.6444444 0.6405229 0.6483660 0.6522876 0.6444444
## [4,] 0.6483660 0.6483660 0.6444444 0.6562092 0.6483660
## [5,] 0.6483660 0.6444444 0.6444444 0.6405229 0.6366013

Color images

print(vesuvio,short=TRUE)
## Image 
##   colorMode    : Color 
##   storage.mode : double 
##   dim          : 851 315 3 
##   frames.total : 3 
##   frames.render: 1

Image stacks

print(nuclei,short=TRUE)
## Image 
##   colorMode    : Grayscale 
##   storage.mode : double 
##   dim          : 510 510 4 
##   frames.total : 4 
##   frames.render: 4

An image histogram shows the distribution of pixel intensities.

hist(greyVesuvio)

range(greyVesuvio)
## [1] 0 1

Displaying

The function display() allows to display images using an interactive JavaScript viewer or R’s built-in graphics capabilities.

display(greyVesuvio)

display(vesuvio)

display(nuclei,all=TRUE)

Writing

writeImage(greyVesuvio,'vesuvio.jpeg',quality = 85)
file.size('vesuvio.jpeg')
## [1] 54781

Image processing

The operations that can be applied to an image can be classified into three categories: point, local and global.

Point operations: output image pixel value at a specific position depends only on the input image pixel value at the same position. Examples: brigthness and contrast adjustment, gamma correction, thresholding.

Local operations: output image pixel value at a specific position depends on the pixel values in the neighborhood of input image pixel value at the same position. Example: convolution.

Global operations: output image pixel value at a specific position depends on all pixel values in the input image.

Brightness adjustment

\(y=x+b\)
display(
   combine(
     vesuvio+0.2,
     vesuvio,
     vesuvio-0.2
   ), all=TRUE
 )

Contrast adjustment

\(y=x*c\)
display(
   combine(
     vesuvio*0.2,
     vesuvio,
     vesuvio*2
   ), all=TRUE
 )

Gamma correction

\(y=x^{\gamma}\)
display(
   combine(
     vesuvio^0.7,
     vesuvio,
     vesuvio^2
   ), all=TRUE
 )

Thresholding

Thresholding is a segmentation method that converts a gray-scale image into a binary image. During this operation all pixels whose intesity values are above a threshold are set to a foreground value, while all the remaining pixels are set to a background value.

\(y = x \le t\) or \(y = x \ge t\)
otsuGreyVesuvio = greyVesuvio < otsu(greyVesuvio)
print(otsuGreyVesuvio)
## Image 
##   colorMode    : Grayscale 
##   storage.mode : logical 
##   dim          : 851 315 
##   frames.total : 1 
##   frames.render: 1 
## 
## imageData(object)[1:5,1:6]
##       [,1]  [,2]  [,3]  [,4]  [,5]  [,6]
## [1,] FALSE FALSE FALSE FALSE FALSE FALSE
## [2,] FALSE FALSE FALSE FALSE FALSE FALSE
## [3,] FALSE FALSE FALSE FALSE FALSE FALSE
## [4,] FALSE FALSE FALSE FALSE FALSE FALSE
## [5,] FALSE FALSE FALSE FALSE FALSE FALSE
display(
  combine(
    greyVesuvio,
    otsuGreyVesuvio
  ),all=TRUE
)

Morphological operations

Morphological operations are usually applied on binary images.

Erosion: the value of the output pixel is the minimum value of all the pixels in the input pixel’s neighborhood.
Dilation: the value of the output pixel is the maximum value of all pixels in the input pixel’s neighborhood.
Opening: is an erosion followed by a dilation.
Closing: is a dilation followed by an erosion.

display(
   combine(
     otsuGreyVesuvio,
     erode(otsuGreyVesuvio),
     dilate(otsuGreyVesuvio),
     opening(otsuGreyVesuvio),
     closing(otsuGreyVesuvio)
   ), all=TRUE
)

fillHull() allows to fill holes in objects.

display(
   combine(
     closing(otsuGreyVesuvio),
     fillHull(closing(otsuGreyVesuvio))
   ), all=TRUE
)

Convolution

Convolution is a mathematical operation on two functions to produce a third function. Convolution can be written in the symbolic form as follows:

\(g(x,y)=(w*f)(x,y)=\sum_{i=-\infty}^{+\infty}\sum_{j=-\infty}^{+\infty}w(i,j)f(x+i,y+j)\)


\(g(x,y)\) is the filtered image, \(f(x,y)\) is the original image, \(w\) is the kernel.

\(I_{1}\) \(I_{2}\) \(I_{3}\) \(I_{4}\) \(I_{5}\) \(I_{6}\) \(I_{7}\) \(I_{8}\) \(I_{9}\)
\(I_{10}\) \(I_{11}\) \(I_{12}\) \(I_{13}\) \(I_{14}\) \(I_{15}\) \(I_{16}\) \(I_{17}\) \(I_{18}\)
\(I_{19}\) \(I_{20}\) \(I_{21}\) \(I_{22}\) \(I_{23}\) \(I_{24}\) \(I_{25}\) \(I_{26}\) \(I_{27}\)
\(I_{28}\) \(I_{29}\) \(I_{30}\) \(I_{31}\) \(I_{32}\) \(I_{33}\) \(I_{34}\) \(I_{35}\) \(I_{36}\)
\(I_{37}\) \(I_{38}\) \(I_{39}\) \(I_{40}\) \(I_{41}\) \(I_{42}\) \(I_{43}\) \(I_{44}\) \(I_{45}\)
\(I_{46}\) \(I_{47}\) \(I_{48}\) \(I_{49}\) \(I_{50}\) \(I_{51}\) \(I_{52}\) \(I_{53}\) \(I_{54}\)
\(I_{55}\) \(I_{56}\) \(I_{57}\) \(I_{58}\) \(I_{59}\) \(I_{60}\) \(I_{61}\) \(I_{62}\) \(I_{63}\)
\(I_{64}\) \(I_{65}\) \(I_{66}\) \(I_{67}\) \(I_{68}\) \(I_{69}\) \(I_{70}\) \(I_{71}\) \(I_{72}\)
\(I_{73}\) \(I_{74}\) \(I_{75}\) \(I_{76}\) \(I_{77}\) \(I_{78}\) \(I_{79}\) \(I_{80}\) \(I_{81}\)
\(K_{1}\) \(K_{2}\) \(K_{3}\)
\(K_{4}\) \(K_{5}\) \(K_{6}\)
\(K_{7}\) \(K_{8}\) \(K_{9}\)

\(O_{12}=I_{2}K_{1}+I_{3}K_{2}+I_{4}K_{3}+I_{11}K_{4}+I_{12}K_{5}+I_{13}K_{6}+I_{20}K_{7}+I_{21}K_{8}+I_{22}K_{9}\)

Build kernel

Low-pass filter

A low-pass filter is used to blur images, remove artefacts. An example of low-pass filter is the gaussian filter. A filter can be generated using makeBrush().

gaussianKern=makeBrush(size = 5, shape = 'gaussian', sigma = 3)
0.0317564 0.0375158 0.0396589 0.0375158 0.0317564
0.0375158 0.0443196 0.0468515 0.0443196 0.0375158
0.0396589 0.0468515 0.0495280 0.0468515 0.0396589
0.0375158 0.0443196 0.0468515 0.0443196 0.0375158
0.0317564 0.0375158 0.0396589 0.0375158 0.0317564
High-pass filter

A high-pass filter is used to detect edges, sharpen images. An example of high-pass filter is the laplacian filter.

values=c(-1,0,-1,0,4,0,-1,0,-1)
laplacianKern=matrix(values,nrow = 3,ncol = 3, byrow = TRUE)
-1 0 -1
0 4 0
-1 0 -1
Laplacian of gaussian filter
values=c(0,0,-1,0,0,0,-1,-2,-1,0,-1,-2,16,-2,-1,0,-1,-2,-1,0,0,0,-1,0,0)
lapgauss5kern=matrix(values,nrow = 5,ncol = 5, byrow = TRUE)
0 0 -1 0 0
0 -1 -2 -1 0
-1 -2 16 -2 -1
0 -1 -2 -1 0
0 0 -1 0 0

Filtering

filter2() allows to filter an image using the 2D FFT (2-Dimensional Fast Fourier Transform) convolution product.

GaussImage=filter2(vesuvio,gaussianKern)
LaplacianImage=filter2(vesuvio,laplacianKern)
LoGImage=filter2(vesuvio,lapgauss5kern)
display (
  combine(
    vesuvio,
    GaussImage,
    LaplacianImage,
    LoGImage
  ), all=TRUE    
)

display (
  combine(
    channel(LaplacianImage,'grey'),
    channel(LoGImage,'grey')
  ), all=TRUE    
)

Example: application in cell biology

“cells.tif” and “nuclei.tif” are fluorescent microscopy images from two channels of HeLa cells, distributed with EBImage.

cells=readImage(system.file("images","cells.tif",package="EBImage"))
nuclei=readImage(system.file("images","nuclei.tif",package="EBImage"))
display(cells,all=TRUE)

display(nuclei,all=TRUE)

rgb = rgbImage(green=1.5*cells,blue=nuclei)
display(rgb,all=TRUE)

Nuclei segmentation

nmaskt = thresh(nuclei, w = 15, h = 15, offset = 0.05)
display(nmaskt,all=TRUE)

nmaskf = fillHull(opening(nmaskt, makeBrush(5, shape='disc')) )
display(nmaskf,all=TRUE)

display( combine(getFrame(nmaskt, 3), getFrame(nmaskf, 3)), all=TRUE )

dmap = distmap(nmaskf)
display(dmap,all=TRUE)

range(dmap)
## [1]  0.00000 14.31782
display(normalize(dmap), frame = 3)

nmask = watershed(dmap, tolerance = 2)
display(nmask,all=TRUE)

display( combine(
  toRGB( getFrame(nuclei, 3) ), 
  colorLabels( getFrame(nmask, 3) )
), all=TRUE )

Cytoplasm segmentation

cmaskt = closing( gblur(cells, 1) > 0.105, makeBrush(5, shape='disc') )
cmask  = propagate(cells, seeds=nmask, mask=cmaskt, lambda = 0.001)
display( combine(
  toRGB( getFrame(cmaskt, 3) ), 
  colorLabels( getFrame(cmask, 3) )
), all=TRUE )

display( paintObjects(nmask,
            paintObjects(cmask, rgb, col = "magenta", thick = TRUE),
         col = "yellow", thick = TRUE), all = TRUE)

st = stackObjects(cmask, rgb)
display(st, all = TRUE)

features=computeFeatures.shape(cmask[,,1], cells[,,1])
head(features)
##   s.area s.perimeter s.radius.mean s.radius.sd s.radius.min s.radius.max
## 1   4216         408      43.52530   14.254634     16.13319     73.48543
## 2   3079         263      34.60167    9.061043     19.14618     53.49633
## 3   3088         211      30.99363    4.630083     16.81697     38.91545
## 4   2377         209      29.37053    8.030861     16.27813     44.85961
## 5   2643         277      34.75310   13.817002     15.03601     62.17141
## 6   2095         181      26.82256    6.948619     16.26868     40.91585

References

Package ‘EBImage’

Introduction to EBImage

Basics of image data and spatial patterns analysis in R